GENERALIDADES

Al igual que el ANOVA y las Correlaciones, la Regresión lineal se basa en una relación LINEAL, pero entre una variable dependiente y una o varias variables independientes.

Busca la relación entre dos variables numéricas, pero tratando de predecir una en función de la otra, es decir, predecir o estimar la variable dependiente a partir de la variable independiente.

El concepto original de Regresión fue introducido por primera vez por Sir Francis Galton a finales del siglo 19, haciendo estudios poblacionales sobre la talla

GENERALIDADES… continuación

Actualmente, el significado de regresión en Estadística, equivale en general al término función en matemáticas, es decir, dado el valor de una variable X, le corresponde un valor de la variable Y.

Esta relación o función lineal, se calcula a partir de un procedimiento llamado rectas de mejor ajuste, con una técnica llamada mímimos cuadrados, originalmente propuesto por Carl Friedrich Gauss. Es decir, la idea es generar una “recta” que mejor se aproxime a los datos de acuerdo con el criterio de mínimo error cuadrático.

Donde:

  • \(\bar{Y}\) = Dependiente
  • X = Independiente
  • a = Coeficiente de Intersección
  • b = Pendiente o Incremento

GENERALIDADES… continuación

¿Como se calculan o se obtienen los coeficientes?

\[b = \frac{\sum_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y})}{\sum_{i=1}^{n}(x_i-\bar{x})^2}\]

\[a = \bar{y}-(b\bar{x})\]

El método de los mínimos cuadrados consiste en hacer un modelo de regresión minimizando la suma de los cuadrados de los errores. Por lo tanto, el criterio de los mínimos cuadrados se basa en minimizar la siguiente expresión:

\[[MIN] \sum_{i=1}^{n}e_i^2=[MIN]\sum_{i=1}^{n}(y_i-\bar{y})^2=[MIN]\sum_{i=1}^{n}(y_i-a-bx_i)^2\]

GENERALIDADES… continuación

El cuadrado de un error es el cuadrado de un residuo, por lo tanto, el cuadrado de un error es igual a la diferencia entre el valor real y el valor ajustado por el modelo de regresión elevado a la dos.

Por eso el criterio de los mínimos cuadrados también se llama criterio del mínimo error cuadrático.

Los modelos de Regresión nos permiten establecer estimaciones o predicciones de determinados valores, con cierta confiabilidad o certeza.

En el análisis de regresión se realizan dos operaciones principales:

  1. Obtener una Ecuación que represente la relación entre las variables (ecuación de regresión) y obtener la línea de regresión dada por esta ecuación. (pueden ser lineales o no)

  2. Estimar una variable llamada variable Dependiente o Respuesta, a partir de otra variable o grupo de variables llamadas variables Independientes o Explicativas.

GENERALIDADES… continuación

Algunos ejemplos de Regresión en Antropología Física:

  • Estimar el peso de una persona a partir de la Talla o Estatura
  • Estimar la estatura de un individuo a partir de la magnitud de una estructura ósea
  • Estimar la masa grasa a partir de pliegues cutáneos
  • Estimar la edad a partir de huesos de la mano y muñeca, rodilla, femur, codo, tobillo, hombro, etc
  • Estimar el sexo a partir de huesos de La pelvis, cráneo y mandíbula, entre otros. (regresion logística)

DIFERENTES MODELOS DE REGRESIÓN

  • Regresión lineal simple –> \(\bar{Y}=a+bX\)
  • Regresión lineal múltiple –> \(\bar{Y}=a+b_1X_1+b_2X_2+...+b_nX_n\)
  • Rgresión no lineal –> Exp: \(\bar{Y}=a*e^{bX}\) , Potencia: \(\bar{Y}=a*X^{b}\) , cuadrático: \(\bar{Y}=a+b_1X+b_2X^2\)
  • Series temporales –> Las series cronológicas son aquellas en que una de las variables en estudio es el tiempo

SUPUESTOS DE LA REGRESIÓN

El modelo de regresión asume ciertos supuestos o condiciones que deben ser tomados en consideración antes de asignarle un valor predictivo

Estos supuestos tienen que ver con la distribución de los residuales Distinguimos tres supuestos principales:

  1. Normalidad de los residuales
  2. Homocedasticidad de los residuales
  3. Independencia de los residuales (o no autocorrelación)

NORMALIDAD DE LOS RESIDUALES

1.- Los residuales deben distribuirse normalmente

Para facilitar la estimación del modelo de regresión es exigible la normalidad de la distribución de los errores.

Se pueden aplicar pruebas de normalidad a la distribución de los residuales, tales como K-S o Shapiro-Wilk. Sin embargo, una alternativa bastante utilizada para este caso es la prueba de Jarque-Bera.

La Prueba de Jarque-Bera parte del hecho de que una distribución normal tiene coeficiente de asimetría igual a 0 y Curtosis igual a 3. Estos dos elementos se miden a partir de los residuos del modelo lineal.

Tambien suelen usarse los gráfico Q-Q Plots o histogramas

En caso de rechazar la normalidad, esta se puede corregir, realizando transformaciones fundamentalmente a la variable dependiente.

Una de las transformaciones más usadas es la Transformación de Box-Cox. Esta es una transformación potencial que corrige la asimetría de una variable, varianzas diferentes o la no linealidad entre variables. En consecuencia, resulta muy útil para transformar una variable y obtener una nueva que siga una distribución normal.

HOMOCEDASTICIDAD DE LOS RESIDUALES

2.- Debe existir homocedasticidad o homogeneidad de las varianzas

El supuesto de homocedasticidad exige que para todo el recorrido de la variable X la varianza del error sea constante, por tanto debe existir homogeneidad de las varianzas de los residuos

Para estos casos se aplica el Test de Breusch-Pagan

Cómo corregir la heterocedasticidad en caso que sea detectada:

  1. Transformar la variable dependiente: Una forma de corregir la heterocedasticidad es transformar la variable dependiente, por ejemplo usando el log.

  2. Redefinir la variable dependiente: Otra forma de corregir la heterocedasticidad es redefinir la variable dependiente. Una forma habitual de hacerlo es utilizar una tasa para la variable dependiente, en lugar del valor bruto.

  3. Utilizar regresión ponderada: Otra forma de corregir la heterocedasticidad es utilizar la regresión ponderada. Este tipo de regresión asigna un peso a cada punto de datos en función de la varianza de su valor ajustado

INDEPENDENCIA DE LOS RESIDUALES

3.- No debe existir autocorrelación entre los residuales

En la regresión lineal y no lineal, se supone que los residuos son independientes unos de otros (no están correlacionados).

Si se viola el supuesto de independencia, algunos resultados de ajustes de modelos pudieran no ser fiables. Por ejemplo, una correlación positiva entre los términos de error tiende a inflar los valores t de los coeficientes, haciendo que los predictores parezcan significativos cuando pudieran no serlo.

Para estos casos se usa el test de autocorrelation Durbin-Watson

En casos que se identifique positivamente la autocorrelación, podemos corregir ese efecto con el metodo de Cochrane-Orcutt paquete “orcutt” en R.

Se ejecuta la función cochrane.orcutt() sobre el modelo de regresión

  • coch1 <- cochrane.orcutt(modelo)
  • coch1

El output muestra como se ajustan los coeficiente y cambia el valor del test Durbin-Watson

SIGNIFICACIÓN DEL MODELO Y R^2

R^2. De los residuos se calcula el Coeficiente de Determinación R2:

Este coeficiente determina Cuánto es la proporción de la varianza de la v. dependiente explicada por las v. independientes.

El R2 va de 0 (el modelo no explica la varianza de la v. dependiente) a 1 (el modelo explica toda la varianza de la v. dependiente).

SIGNIFICACIÓN DEL MODELO

El modelo de regresión también presenta un valor F y un p-value, que representa la significación del modelo completo.

EJEMPLO

Usaremos la variable Peso en Kilogramos (wgt) como independiente y la variable talla (ht) en cms de la base base_enca2014_200.xlsx

Primero calculamos la correlación para tener idea del nivel de asociación entre las variables

Observamos que la correlación entre ambas es de 0.445 y a pesar que podemos considerarla moderada, es significativa porque su p-valor es menor a 0.05

EJEMPLO… continuación

Procedemos a calcular los coeficientes de la regresion lineal

\[b=\frac{11616.4}{41011.3}=0.283 \hspace{2em} a=158.4-(0.283*72.2)=138.210\]

\[Ecuación \hspace{0.5em} de \hspace{0.5em} Regresión: \hspace{0.5em}y = 158.4+(0.283*x) \]

EJEMPLO… continuación

Revisamos la salida de R

EJEMPLO… continuación

A partir de la Suma Cuadrados de los residuos: \(\sum(y-y_e)^2\) y Total: \(\sum y^2-(\frac{(\sum y)^2}{n})\) , podemos construir la tabla de ANOVA del modelo

Observamos que el valor F, es el mismo que arroja la salida de R, así como su p-valor

EJEMPLO… continuación

Verificamos los supuestos

EJEMPLO… continuación

Revisamos gráficamente

EJEMPLO… continuación

Otra forma de ver los residuos

EJEMPLO… continuación

Dado que se ha detectado la heterocedasticidad, intentaremos realizar el ajuste ponderado.

Específicamente realizaremos los mínimos cuadrados ponderados definiendo los pesos de tal manera que se dé más peso a las observaciones con menor varianza.

Es decir, realizaremos una nueva regresión donde se ponderan las observaciónes en función de las varianzas.

Vemos como se ha corregido la heterocedasticidad al correr nuevamente el test de Breusch-Pagan en el nuevo modelo de regresión ponderado

EJEMPLO… continuación

Revisamos el nuevo modelo

A partir de la salida, podemos ver que la estimación del coeficiente para la variable predictora wgt cambió un poco y el ajuste general del modelo mejoró.

El modelo de mínimos cuadrados ponderados tiene un error estándar residual de 1,19 en comparación con 8.502 en el original simple modelo de regresión lineal.

Esto indica que los valores predichos producidos por el modelo de mínimos cuadrados ponderados están mucho más cerca de las observaciones reales en comparación con los valores predichos producidos por el modelo de regresión lineal simple.

El modelo de mínimos cuadrados ponderados también tiene un R-cuadrado de 0.19 en comparación con 0.18 en el modelo de regresión lineal simple original.

EJEMPLO… continuación

Finalmente realizamos el gráfico de dispersión

PRACTICA

PRACTICA… continuación

Usaremos una muestra aleatoria de tamaño 200 de la Encuesta nacional de Consumo Alimentario del 2014, ya trabajadas en clases anteriores (base_enca2014_200.xlsx)

Usaremos la variable Peso en Kilogramos (wgt) como independiente y la variable talla (ht) en cms

Primero calculamos la correlación para tener idea del nivel de asociación entre las variables

#library(readxl)
#datos <- read_excel("data/base_enca2014_200.xlsx")
#cor1 <- cor.test(datos$ht, datos$wgt, method = "spearman", exact = FALSE)
cor1
    Spearman's rank correlation rho

data:  datos$ht and datos$wgt
S = 739489, p-value = 3.893e-11
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.4453693 

Observamos que la correlación entre ambas es de 0.445 y a pesar que podemos considerarla moderada, es significativa porque su p-valor es menor a 0.05

PRACTICA… continuación

Procedemos a calcular la regresion lineal: reglineal1 <- lm(ht~wgt, data=datos)

#reglineal1 <- lm(ht~wgt, data=datos)
summary(reglineal1) 
Call:
lm(formula = ht ~ wgt, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
-18.791  -6.289  -1.368   5.846  22.130 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 138.21023    3.04934  45.325  < 2e-16 ***
wgt           0.28325    0.04198   6.747 1.63e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.502 on 198 degrees of freedom
Multiple R-squared:  0.1869,    Adjusted R-squared:  0.1828 
F-statistic: 45.52 on 1 and 198 DF,  p-value: 1.626e-10

PRACTICA… continuación

Podemos “adornar” la salida con el paquete Flextable, tambien hay que usar el dplyr

#library(dplyr)
#library(flextable)
reglineal1 %>% as_flextable() %>% bg(bg="#2874A6", part = "header") %>% 
  bold(part = "header") %>% color(color = "white", part = "header") %>% autofit()

PRACTICA… continuación

Verificamos los supuestos

#library(fBasics)     #Para el test jarqueberaTest()
#library(lmtest)      #Para los test dwtest() y bptest()
jarqueberaTest(reglineal1$residuals)    # test Jarque--Bera de normalidad
Title:
 Jarque - Bera Normalality Test

Test Results:
  STATISTIC:
    X-squared: 4.899
  P VALUE:
    Asymptotic p Value: 0.08634 

Description:
 Mon Oct 23 21:55:10 2023 by user: eluis
dwtest(reglineal1)                      # test de autocorrelation Durbin-Watson
    Durbin-Watson test

data:  reglineal1
DW = 2.1118, p-value = 0.7872
alternative hypothesis: true autocorrelation is greater than 0
bptest(reglineal1)                      # test de heterocedasticidad de Breusch-Pagan
    studentized Breusch-Pagan test

data:  reglineal1
BP = 12.666, df = 1, p-value = 0.0003724

PRACTICA… continuación

par(mfrow=c(2,2))
plot(reglineal1)

par(mfrow=c(1,1))

PRACTICA… continuación

Otra forma de ver los residuos

plot(reglineal1$residuals, type="l", ylab="Residuos", xlab="Observación", 
     col="blue", lwd=3, main="Gráfico de dispersión de los residuos")
abline(h=0, col="black", lwd=2)

PRACTICA… continuación

Dado que se ha detectado la heterocedasticidad, intentaremos realizar el ajuste ponderado.

#Primero definimos la ponderación y la guardamos en un objeto llamado "w"
#w <- 1/lm(abs(reglineal1$residuals)~reglineal1$fitted.values)$fitted.values^2  

#Usamos esa ponderación un un nuevo modelo de regresión y lo incorporamos en "weights"
#reglineal1a <- lm(ht~wgt, data=datos, weights = w)
bptest(reglineal1a)   #Probamos el modelo con el test de Breusch-Pagan
    studentized Breusch-Pagan test

data:  reglineal1a
BP = 1.5198, df = 1, p-value = 0.2176

Vemos como se ha corregido la heterocedasticidad al correr nuevamente el test de Breusch-Pagan en el nuevo modelo de regresión ponderado. Es decir, anteriormente el test arrojaba un p-valor de 0.0003724 y ahora con el nuevo modelo ponderado el test tiene un p-valor de 0.2176. Por tanto, al ser mayor de 0.05 se acepta la homocedasticidad

PRACTICA… continuación

Revisamos el nuevo modelo

summary(reglineal1a)
Call:
lm(formula = ht ~ wgt, data = datos, weights = w)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-2.8328 -0.8910 -0.1895  0.8676  2.8397 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 136.78661    2.92049  46.837  < 2e-16 ***
wgt           0.30351    0.04295   7.067 2.64e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.19 on 198 degrees of freedom
Multiple R-squared:  0.2014,    Adjusted R-squared:  0.1974 
F-statistic: 49.94 on 1 and 198 DF,  p-value: 2.643e-11

PRACTICA… continuación

Vamos a adornar los resultados

reglineal1a %>% as_flextable() %>% bg(bg="#2874A6", part = "header") %>% bold(part = "header") %>% color(color = "white", part = "header") %>% autofit()

A partir de la salida, podemos ver que la estimación del coeficiente para la variable predictora wgt cambió un poco y el ajuste general del modelo mejoró.

PRACTICA… continuación

Finalmente realizamos el gráfico de dispersión

#library(ggplot2)
ggplot(datos, aes(x =wgt , y = ht)) + geom_point() + 
  stat_smooth(method = "lm", col = "red", se = FALSE) + theme_bw() + 
  geom_text(x = 50, y = 180, label = paste("r: ", round(cor1$estimate, 3)), col="blue"  )

OTRO EJEMPLO… continuación

Realizaremos una regresion lineal múltiple, donde consideraremos como dependiente a la variable colest_mg y como independientes a cho_g_dia,prot_g_dia y lip_g_dia.

Primero hacemos una matriz de correlación entre las variables involucradas

cor2 <- cor(datos[, c("colest_mg", "cho_g_dia","prot_g_dia", "lip_g_dia")], method = "spearman")
library(corrplot)
corrplot(cor2, method ="number", tl.col = "black", tl.srt = 20,
         col=colorRampPalette(c("blue","lightblue","red"))(100))

OTRO EJEMPLO… continuación

Hay revisar las variables

summary(datos[ , c("colest_mg", "cho_g_dia","prot_g_dia", "lip_g_dia")])
   colest_mg        cho_g_dia         prot_g_dia       lip_g_dia     
 Min.   :  0.00   Min.   :  55.61   Min.   :  7.42   Min.   :  0.94  
 1st Qu.: 77.93   1st Qu.: 162.37   1st Qu.: 35.91   1st Qu.: 31.65  
 Median :138.84   Median : 233.22   Median : 48.10   Median : 45.06  
 Mean   :195.10   Mean   : 253.33   Mean   : 55.80   Mean   : 51.69  
 3rd Qu.:279.86   3rd Qu.: 314.69   3rd Qu.: 67.66   3rd Qu.: 63.97  
 Max.   :858.56   Max.   :1250.94   Max.   :227.03   Max.   :252.18  

Vemos que en colest_mg, hay un valor que es igual a 0. por tanto procedemos a eliminarlo de la base, porque eso nos puede generar una distorsión en los resultados

datos <- datos[datos$colest_mg>0,]
#Revisamos nuevamente
summary(datos$colest_mg)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.47   79.31  139.50  196.08  281.21  858.56 

Observamos que lo hemos corregido

OTRO EJEMPLO… continuación

Realizamos entonces la regresión

reglineal2 <- lm(colest_mg~cho_g_dia+prot_g_dia+lip_g_dia, data=datos)
#Revisamos los supuestos
jarqueberaTest(reglineal2$residuals) # test Jarque--Bera de normalidad
Title:
 Jarque - Bera Normalality Test

Test Results:
  STATISTIC:
    X-squared: 153.9961
  P VALUE:
    Asymptotic p Value: < 2.2e-16 

Description:
 Mon Oct 23 21:55:11 2023 by user: eluis
dwtest(reglineal2) # test de autocorrelation Durbin-Watson
    Durbin-Watson test

data:  reglineal2
DW = 2.0813, p-value = 0.7141
alternative hypothesis: true autocorrelation is greater than 0
bptest(reglineal2) # test de heterocedasticidad de Breusch-Pagan
    studentized Breusch-Pagan test

data:  reglineal2
BP = 4.2716, df = 3, p-value = 0.2336

OTRO EJEMPLO… continuación

Como la distribución de la normalidad en los residuos se ha rechado, intentaremos realizar una transformación usando el método de Box-Cox.

Vamos a crear una nueva variable Y (dependiente) elevada a un exponente llamdo lambda, generado por el método de Box-Cox.

#install.packages("MASS")
library(MASS)
bc <- boxcox(lm(datos$colest_mg ~ 1))

OTRO EJEMPLO… continuación

Observamos que el valor de lambda está entre 0 y 1, pero para obtener el valor exacto, usamos el siguiente código:

lambda <- bc$x[which.max(bc$y)]
lambda
[1] 0.3030303

Luego lo usamos para transformar la variable dependiente y lo guardamos en nuestra base

datos$colest_mgt <- datos$colest_mg^lambda

Finalmente repetimos la regresión usando la nueva variable transformada

reglineal2a <- lm(colest_mgt~cho_g_dia+prot_g_dia+lip_g_dia, data=datos)
#Verificamos la normalidad
jarqueberaTest(reglineal2a$residuals)
Title:
 Jarque - Bera Normalality Test

Test Results:
  STATISTIC:
    X-squared: 6.1445
  P VALUE:
    Asymptotic p Value: 0.04632 

Description:
 Mon Oct 23 21:55:12 2023 by user: eluis

OTRO EJEMPLO… continuación

Observamos que, aunque lo ajustamos, no supera el valor 0.05, por tanto se mantiene la no normalidad.

Aún así vamos a interpretar la salida de la regresión

reglineal2a %>% as_flextable() %>% bg(bg="#2874A6", part = "header") %>% bold(part = "header") %>% color(color = "white", part = "header") %>% autofit()

Observamos que la variable cho_g_dia, no contribuye con el modelo de estimación. El R2 es 0.42 quiere decir, que la variabilidad del colest_mg es explicada en un 42% por cho_g_dia, prot_g_dia, lip_g_dia. El valor F es significativo, es decir, el modelo explica la estimación. La Ecuación predictiva es: \(Y = 3.021+(-0.001*cho_g_dia)+(0.017*prot_g_dia)+(0.017*lip_g_dia)\)